Phase Diagram for the Winfree Model of Coupled Nonlinear Oscillators 
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In 1967 Winfree proposed a mean-field model for the spontaneous synchronization of chorusing 
crickets, flashing fireflies, circadian pacemaker cells, or other large populations of biological oscil- 
lators. Here we give the first bifurcation analysis of the model, for a tractable special case. The 
system displays rich collective dynamics as a function of the coupling strength and the spread of 
natural frequencies. Besides incoherence, frequency locking, and oscillator death, there exist novel 
hybrid solutions that combine two or more of these states. We present the phase diagram and derive 
several of the stability boundaries analytically. 



PACS numbers: 05.45.Xt, 87.10.+e 

The collective behavior of limit-cycle oscillators was 
first investigated by Winfree Q|. Using a mean- field 
model of coupled phase oscillators with distributed natu- 
ral frequencies, he discovered that collective synchroniza- 
tion is a threshold phenomenon, the temporal analogue of 
a phase transition. Specifically, when the strength of the 
coupling exceeds a critical value, some oscillators sponta- 
neously synchronize to a common frequency, overcoming 
the disorder in their natural frequencies. The model was 
subsequently refined by Kuramoto |i| and others HQ, 

ncu- 
bub- 



with applications to Josephson junction arrays 
trino flavor oscillations j6|, Brownian ratchets ^ 
bly fluids ^ , semiconductor laser arrays ]9| and Landau 
damping of plasmas jl0| . 

Despite all the activity that Winfree's work ultimately 
provoked, surprisingly little is known about the dynamics 
of his original model. In this Letter, we explore a special 
case of the model for which several analytical results can 
be obtained. This version of the model is also related 
to recent work on pulse-coupled oscillators (where the 
oscillators interact by firing sudden impulses), a case of 
interest in neurobiology ftM- However, we are motivated 
here by a dynamical systems perspective: the goal is to 
understand the collective behavior and bifurcations of 
the model as a function of two parameters, the coupling 
strength and the spread of natural frequencies. 

In the limit of weak coupling and nearly identical fre- 
quencies, our model reduces to the Kuramoto model, 
whose behavior is well understood it displays 

locked, partially locked, or incoherent states, depending 
on the choice of parameters. Away from this familiar 
regime, we discover novel hybrid states corresponding to 
various mixtures of locking, incoherence, and oscillator 
death (a cessation of oscillation caused by excessively 
strong coupling JT^j). 

The Winfree model is 
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for i = 1, . . . , N, where N 3> 1. Here 0j(f) is the phase of 
the ith oscillator at time t, k > is the coupling strength, 
and the frequencies u>i are drawn from a symmetric, uni- 



modal density g(w). We assume that the mean of g(ui) 
equals 1, by a suitable rescaling of time, and that its 
width is characterized by a parameter 7. The coupling 
in (|l|) has the following interpretation. The jth oscilla- 
tor makes its presence felt through an influence function 
P(0j); in turn, the ith. oscillator responds to the mean 
field (the average influence of the whole population) ac- 
cording to a sensitivity function R(9i). 

From now on, we consider the special case where 



P{9) = 1 + cos0, R(9) = -sml 



(2) 



Note that this P(6) is a smooth but pulse-like func- 
tion. (At the end of this paper we consider a much 
more sharply peaked P(9); then 9 = represents the 
phase when the oscillator suddenly fires.) The func- 
tional form of R{9) is inspired by the qualitative shape 
of the phase-response curve of many biological oscilla- 
tors |l2|,[l3|]. With these choices, Eq. (|l]) becomes a simple 
model for a population of pulse-coupled biological oscil- 
lators, such as crickets 14 1, fireflies ]l5[ ], or heart pace- 
maker cells [|l6| . The unusual aspect is that the coupling 
is through the phase-response curve JT^ ] ; thus an oscilla- 
tor can be either advanced or delayed by a pulse from an- 
other oscillator, depending on its phase when it receives 
the stimulus. This differs from the strictly excitatory or 
inhibitory coupling often used in integrate-and-fire mod- 
els of neural oscillators. 

We begin by describing our numerical simulations. 
Equation (Q) was integrated numerically using N = 800 
oscillators. The frequencies a;, where chosen to be evenly 
spaced in the interval / = [1 — 7, 1 +7], corresponding to 
a uniform density g(u>) = I/27 for u> €. I, and g(u>) = 
otherwise. The long-term behavior of the system was 
always found to be independent of the initial conditions. 

To compare the long-term dynamics of individual 
oscillators, we compute the average frequency (also 
known as the rotation number) of each oscillator, 
Pi = limt—joo 9i(t)/t. In our simulations, the limit was 
typically well-approximated by integrating the system for 
500 time units, starting from a random initial condition, 
although longer runs were sometimes necessary. The ro- 
tation numbers provide a convenient measure of synchro- 
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FIG. 1. Collective states, as indicated by rotation numbers 
for k = 0.65. Equation with P(9), R(0) as in (0) was 
integrated for 500 time units starting from a random initial 
condition, (a) 7 = 0.1: locking, (b) 7 = 0.205: partial 
locking, (c) 7 = 0.3: incoherence, (d) 7 = 0.6: partial death. 

nization: two or more oscillators are frequency locked if 
they have the same rotation number. 

Figure 1 plots pi as a function of u>i for increasing val- 
ues of 7, at a fixed k. For small 7, all the oscillators 
are locked [Fig. 1(a)]. They can be visualized as a pack 
of particles rotating at the same average rate around the 
unit circle, where Oi(t) denotes the angular position of os- 
cillator i. As 7 is increased past a critical threshold, the 
coupling is no longer sufficient to keep all the oscillators 
mutually entrained. Just above threshold, the system 
stays partially locked, with the fastest oscillators peeling 
away from the pack, but drifting incoherently relative to 
one another [Fig. 1(b)]. With further increases in 7, suc- 
cessively more oscillators peel away until eventually the 
entire population is incoherent [Fig. 1(c)]. For sufficiently 
large 7, the system converges to a state of partial death 
in which the slowest oscillators stop moving altogether, 
while the faster ones remain incoherent [Fig. 1(d)]. 

Partial locking and partial death are hybrid states; 
their rotation number plots [Figs. 1(b) and (d)] contain 
two distinct branches that correspond to qualitatively 
different dynamics. We have also observed hybrid states 
with three and four branches [Fig. 2]. These more com- 
plicated states should all be regarded as variants of par- 
tial locking, since there is at least one branch of locked 
oscillators in each case. (Note that all these states are 
near each other in parameter space.) We label them ac- 
cording to their branches, as follows. Locked- slipping- 
locked [Fig. 2(a)]: There are two separate plateaus of 
locked oscillators, at close to 2:1 frequency ratio in the 
example shown, separated by a branch of slipping os- 
cillators. A slipping oscillator typically co- rotates with 
a locked group for several periods, then slips away for a 
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FIG. 2. Partially-locked hybrid states. (a) 7 = 0.19, 
K = 0.78: locked-slipping-locked. (b) 7 = 0.21, k = 0.76: 
slipping- locked- incoherent, (c) 7 = 0.205, k = 0.79: quiver- 
ing-slipping- locked- incoherent . 

few cycles before eventually rejoining the same group and 
repeating the pattern. Oscillators slip more or less fre- 
quently depending on their native frequency u>i . Slipping- 
locked-incoherent [Fig. 2(b)]: There is a central group of 
locked oscillators, flanked by slower ones that slip and 
faster ones that drift monotonically. Quivering- slipping- 
locked-incoherent [Fig. 2(c)]: This state exists near par- 
tial death. It is similar to the state in Fig. 2(b), but with 
an added mode of behavior: the slowest oscillators quiver 
about their former death phases. An oscillator that quiv- 
ers has zero rotation number — it remains trapped in the 
neighborhood of a single phase for all time — and hence is 
effectively dead, though not completely motionless. Fig- 
ure 3 summarizes the system's long-term behavior as a 
function of k and 7. 

We now outline our analytical calculations of the 
boundaries of the regions corresponding to incoherence, 
partial death, and death. Following the standard ap- 
proach used for the Kuramoto model we rewrite 
the dynamics in the infinite- A limit. By analogy with 
the continuum limit in fluid mechanics, we view the os- 
cillators as particles moving around the unit circle. For 
each frequency to, let p (9, t, u) denote the density of os- 
cillators at phase 6 at time t, and let v(6,t,ui) denote 
the local velocity field. Then p satisfies the continuity 
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FIG. 3. Phase diagram for Eqs. (^) and (^), assuming a uni- 
form distribution of natural frequencies on [1 — 7, 1 +7]. The 
boundary between locking and partial locking is determined 
numerically; all other boundaries are determined analytically. 
All of the partially-locked hybrid states are lumped together 
in one region, for simplicity. 



equation dp/dt — —d(pv) /d6, expressing conservation 
of oscillators of frequency to. The velocity v is deter- 
mined by applying the law of large numbers to Eq. ([!]) . 
The sum over all oscillators in (Q) is replaced by an in- 
tegral as N — > 00, yielding v{9 1 t,uj) = to — a(t)sm9, 
where 



a(t) 



2tt 
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(1 + cas6)p(6,t,u)g(to)(kod6. (3) 
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Additionally, we demand that p be non-negative, 
27r-periodic in 9, and we impose the normalization 
f^p(9,t,w)d9 = 1 for all t,to. 

The key to the analysis is recognizing that incoher- 
ence, partial death, and death correspond to stationary 
densities po(9, to) in the infinite- N limit. Hence, one may 
solve for all three states by seeking fixed points of the 
continuity equation. These satisfy pqVq = C(oj), with 
C(to) determined by normalization. Depending on its 
natural frequency, an oscillator's steady-state behavior 
falls into one of two categories, (i) to < a: In this case, 
v = to — a - sin 6* = 0. The oscillators of a given fre- 
quency to < a all remain stuck at a single phase 6* (to), 
defined implicitly by sin 9* = to /a. Their density is 
Po(9,to) = 5(9 — 9*(to)). Such oscillators are motion- 
less, or dead, (ii) to > a: These oscillators rotate non- 
uniformly around the circle, hesitating near 9 = n/2 and 
accelerating near 9 = 3ir/2, as dictated by their veloc- 
ity field vq = to — a sin 9 > 0. The stationary density is 
inversely proportional to the velocity: 
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uj — a sin ( 



(4) 



Normalization then implies C(to) — \fto 2 — a 2 / (2ir). 
From these two basic scenarios, we are able to calculate 



the incoherence, partial death and death boundaries as 
follows. 

Incoherence exists provided to > a for all to; then all 
oscillators belong to category (ii) above. The bound- 
ary separating incoherence and partial death occurs when 



where to n 



1 — 7. The first oscillators 



to die are the slowest ones. To solve for a, we substitute 
into Eq. (||); this yields a = k. Thus, partial death 
bifurcates from incoherence along the straight line 



K= 1 
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assuming k is not so large that the system is in the death 
region. Remarkably, this result holds for any frequency 
distribution g(to), whether symmetric or not. 

The stability of the incoherent state is determined by 
linearizing the continuity equation about the incoherent 
density (0). The resulting linear operator has a continu- 
ous spectrum that is pure imaginary and a discrete spec- 
trum that is governed by the equation |lq] 
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A 2 + UJ 2 - K 2 



■ g(oj)doj. 
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From (^), it is clear there are no eigenvalues A with 
Re(A) < 0; if there were, the right-hand side would have 
negative real part, contradicting the assumption n > 0. 
Thus, incoherence is either unstable or neutrally stable. 
Numerics indicate that the boundary between incoher- 
ence and partial locking corresponds to a Hopf bifurca- 
tion. To obtain the boundary, we solve Eq. @ for A 
perturbatively, assuming 7 <C 1, and then take the limit 
Re(A) — > + . The result for a uniform g(to) is 
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+ 0( 7 7 ). (7) 



In the limit 7 — > 0, Eq. (0) reduces to k = 87/71", which 
is the critical coupling threshold for the averaged system 
(the Kuramoto model, with coupling K = k/2). Hence, 
Eq. (|7]) can be viewed as an extension of the classical 
threshold condition [QJ2) into the non-averaged regime of 
stronger coupling and frequency disorder. 

Finally, to calculate the death boundary we use a self- 
consistency argument familiar from mean- field theory pi . 
In the death state, each oscillator comes to rest at 6* (to), 
defined by the zero velocity condition sin#* = to /a. 
This requires a > uj max , where uj max = 1+7. Each 
phase 9* (to) depends on cr, which in turn depends on 
all phases via Eq. (Q); therefore, a must be determined 
self-consistently. For each to, there are two possible 
roots 9*(to): one in [0,7r/2], the other in [7r/2,7r]. How- 
ever, the unique stable fixed point of Eq. ([[]) satisfies 
< 9* (to) < it/2 for all to. Substituting the correspond- 
ing density po(9, to) = 8(9 — 9* (to)) into Eq. (||) gives the 
self-consistency condition 
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FIG. 4. Phase diagram for (|lj) with P(6) = a„(l + cos6») n , 
n — 10 and a uniform frequency distribution on [1 — 7, 1 + 7]. 
Death boundary determined analytically; all other boundaries 
determined numerically. 



For death to exist, there must be a root a > 1 + 7 of 
Eq. (g). The boundary between death and partial death 
corresponds to an endpoint bifurcation, and is found by 
setting <r = 1 + 7 in For a uniform g(u>), this yields 
the exact expression 
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where A = (1 — 7)/(l+7). The remaining portion of the 
boundary, separating death from full and partial locking, 
corresponds to a saddle-node bifurcation. We obtain it by 
solving Eq. (||) numerically, together with the tangency 
condition 1/k = F'(a), where F(a) is the right-hand side 
Of Eq. (§). 

To check the robustness of the phase diagram, we re- 
placed P(0) in (|^) with a family of influence functions 
Pn{9) = c n (l + cos#) n , n > 1 which becomes more 
and more sharply peaked as n increases. (The normal- 
ization coefficients a n are determined by requiring P n (0) 
to have integral equal to 2n over one cycle. Note that 
P n {6) -> 2ttS(9) as n -> 00.) We find that all of the 
phenomena observed for the model studied in this paper 
(n = 1) persist as we increase n; the only difference is 
that the boundaries become slightly distorted [Fig. 4]. 

The mean-field model ([!]) is one of the simplest pos- 
sible models of pulse-coupled oscillators. More realistic 
models would include such features as spatial coupling, 
time delay, dynamical synapses, refractory period, non- 
sinusoidal influence and sensitivity functions, etc. It re- 
mains to be seen whether such models would also exhibit 
the hybrid states found here. In any case, we now know 
that even the most idealized version of the Winfree model 
displays a fascinating wealth of dynamics that, curiously, 
escaped notice for over thirty years. 
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